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We investigate the properties in finite magnetic field of an extended anisotropic XXZ spin-1/2 
model on the Kagome lattice, originally introduced by Balents, Fisher, and Girvin [Phys. Rev. 

B, 65, 224412 (2002)]. The magnetization curve displays plateaus at magnetization m = 1/6 and 
1/3 when the anisotropy is large. Using low-energy effective constrained models (quantum loop and 
quantum dimer models), we discuss the nature of the plateau phases, found to be crystals that break 
discrete rotation and/or translation symmetries. Large-scale quantum Monte-Carlo simulations 
were carried out in particular for the m = 1/6 plateau. We first map out the phase diagram of the 
effective quantum loop model with an additional loop-loop interaction to find stripe order around 
the point relevant for the original model as well as a topological Z 2 spin liquid. The existence 
of a stripe crystalline phase is further evidenced by measuring both standard structure factor and 
entanglement entropy of the original microscopic model. 
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I. INTRODUCTION 

The study of the disruption of the ordering tendency 
in low-dimensional antiferromagnets at zero tempera¬ 
ture and the subsequent emergence of unconventional 
phases is an active topic of research, on both theoretical 
or experimental fronts. Frustrating interactions, quan¬ 
tum fluctuations or imposition of a magnetic field are 
among the main ingredients to destabilize the antiferro¬ 
magnetic long-range order [Ij. Considerable efforts have 
been devoted to studying and characterizing the result¬ 
ing phases of matter, which may assume another type of 
non-magnetic ordering (such as crystalline states of sin¬ 
glet valence bonds), or may not exhibit any kind of local 
order at all {quantum spin liquids [2]). 

Two-dimensional non-bipartite lattice with antiferro¬ 
magnetic interactions between spin-1/2 moments are 
amongst the most studied, as antiferromagnetic order is 
most severely challenged by quantum fluctuations there. 
The recent numerical findings of topological spin liquid 
phases in the realistic SU(2) Kagome-lattice antiferro¬ 
magnet [3-5], or its XXZ version [G] provide examples of 
stabilization of such exotic phases of matter. Previously, 
Z 2 topological spin liquids were found in toy models, such 
as constrained quantum dimer models (QDM) [7] on the 
frustrated triangular [8] or Kagome [9] lattices. These 
simplified models advantageously allow for a better an¬ 
alytical handle and understanding of the physics of spin 
liquids, but their connection to SU(2) microscopic Hamil¬ 
tonians is not always direct. 

Another route to realizing exotic phases of matter 
is to start from a model with highly degenerate man¬ 
ifold of ground states {ice manifold) obeying local con¬ 
straints and then derive an effective Hamiltonian describ¬ 
ing the emergent degrees of freedom in this manifold [10- 
12]. Following this strategy, Balents, Fisher and Girvin 


(BFG) introduced a spin-1/2 XXZ model with extended 
interactions on the Kagome lattice [10]. Using a general¬ 
ized QDM derived for low energies, these authors argued 
that this system hosted a topological gapped Z 2 spin liq¬ 
uid phase [10], which was further confirmed numerically 
in subsequent works [13-15]. One of the hallmarks of a 
Z 2 liquid phase is the existence of a topological correction 
to entanglement entropy, which was shown to be present 
in the BFG model [16]. 

In considering the nature of the (gapped) phases re¬ 
sulting from the destruction of antiferromagnetic or¬ 
der, one important guiding principle is provided by a 
general statement about the relation between the pres¬ 
ence/absence of an excitation gap and commensurabil- 
ity. For instance, it can be rigorously shown [17, 18] that 
a given spin system (in the absence of magnetic field) 
is gapless when the value of total spin per unit cell is 
an half-odd-integer (this is the case for, e.g., the spin- 
1/2 Heisenberg model on the Kagome lattice) and the 
ground state is non-degenerate. That is, if we have a fea¬ 
tureless liquid-like ground state with an excitation gap, 
the state necessarily has a hidden degeneracy (probably 
of topological origin). A simpler version of this argument 
[19] applies to any spin-^ systems with U(l) symme¬ 
try and is of direct relevance to magnetization plateaux 
(see [1] for a review): a unique ground state with a gap 
is possible only when the number of spins q within a 
lattice unit cell and the magnetization m per site sat¬ 
isfy q{S — m) € Z. Recently, a field-theoretical mean¬ 
ing was given to this relation [20] and it was predicted 
that when the q{S — m) is a simple (non-integral) ra¬ 
tional number, we may have gapped featureless liquid 
phases {i.e., spin-liquid plateaus) as well as more con¬ 
ventional plateaus with magnetic superstructures. The 
search for microscopic models that potentially host these 
spin-liquid plateaus formed in high magnetic fields is 
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quite intriguing in the light of both the field control of 
spin liquids [21] and a recent theoretical report [22]. 

In this paper, we will investigate the ground-state 
phases of the BFG model in the presence of a magnetic 
field, which, as we will see, exhibits large magnetization 
plateaus at magnetization per site m = 1/6 and 1/3 as 
well as the one at m = 0 (the gapped spin liquid reported 
in Ref. [10]). On these plateaus, the low-energy physics 
is also well captured by effective constrained models on 
a triangular lattice: a generalized QDM for the m = 0 
plateau [10], the usual QDM for m = 1/3, and a quan¬ 
tum loop model (QLM) for m = 1/6. While much is 
known already about the nature of the to = 1/3 plateau 
thanks to numerous extensive studies [8, 23-29] done for 
the QDM on the triangular lattice, the phase diagram of 
the QLM, which is of direct relevance to the to = 1/6 
plateau, is not known to our best knowledge. We will 
investigate it in details using simple analytical consider¬ 
ations and quantum Monte Carlo (QMC) simulations. 

The paper is organized as follows. In Sec. II, we present 
the spin 1/2 model considered in this study and provide 
the low-energy mapping to constrained models valid for 
the magnetization plateaus at to = 1/6 (QLM) and 1/3 
(QDM). We also review some key results obtained pre¬ 
viously when no magnetic field is present. Sec. Ill is 
devoted to a complete mapping of the phase diagram of 
the QLM on the triangular lattice. For the QLM with 
only the kinetic term, that is relevant to the original spin 
model, we find that the ground state of the QLM displays 
crystalline behavior with a spontaneous breaking of six¬ 
fold rotation symmetry, but not of translation symme¬ 
try. We will then turn, in Sec. IV, to large-scale numeri¬ 
cal simulations of the original microscopic spin model at 
TO = 1/6 using direct QMC simulations. Computation 
of the diagonal spin structure factor confirms the exis¬ 
tence of a stripe ordering of down spins with a moderately 
large correlation length. The presence of this crystalline 
phase is further corroborated using Renyi entanglement 
entropy. Conclusions are given in Sec. V. 


II. EXTENDED XXZ MODEL IN A MAGNETIC 
FIELD ON THE KAGOME LATTICE 

A. Microscopic model 

In the search for simple microscopic models (e.g. 
with two-spin interactions) hosting quantum spin liquid 
phases, BFG introduced a spin-1/2 XXZ model on the 
Kagome lattice [10] with the aim to reproduce, in a cer¬ 
tain regime of parameters, the physics of a (generalized) 
quantum dimer model, known for hosting a Z 2 topologi¬ 
cal spin liquid phase. We reproduce here their construc¬ 
tion for completeness and obtain the effective Hamiltoni¬ 
ans for the plateaus. The Hamiltonian is defined on the 
hexagonal plaquettes of the Kagome lattice (see Fig. 1) 



FIG. 1. (Color online) Kagome lattice with the interactions 
between the first (black), second (red) and third (blue) neigh¬ 
bours within the same hexagonal plaquettes, as included in 
the model (1). Typical examples of the third-nearest-neighbor 
interaction between different hexagons [that are not included 
in the Hamiltonian (1)] are shown by green dashed lines. Oth¬ 
ers are obtained by space-group operations. 


and reads: 


H — Hq + Hxy, 

= X X sts--kj2st , , 

O L (bi>60 ieO J ^ 

H,y=Y,J.y Y. {StS-+S-S+), 

O <bi>60 

where is the spin-1/2 operator. Hq is the diagonal 
part of the Hamiltonian made up of the Ising interac¬ 
tions between first, second and third neighbors within 
each plaquette (see Fig. 1 for all the corresponding links) 
as well as the Zeeman term for a magnetic field h along 
the z-axis. Note that the third-neighbor interactions be¬ 
tween sites belonging to different hexagons, shown by the 
green dashed lines in Fig. 1, are not included (see Sec. HB 
for the effect of these neglected interactions). The sec¬ 
ond part H,cy contains the spin-flip terms between sites 
within the same hexagon. Originally, BFG considered all 
such terms on each plaquette, but, as will be clear be¬ 
low, setting J^y = 0 for the second and third neighbors 
does not alter the physical behavior [see Eq. (4)]. Com¬ 
pared to the original Hamiltonian defined in Ref. 10, we 
added the magnetic field term and the coefficient of the 
Jz terms is modified by a factor of 1/2. The model can 
also be reformulated in terms of hardcore bosons, with 
the boson density playing the role of the magnetization. 
In this language, the plateaus at to = 1/6 and to = 1/3 
that we will investigate correspond to 1 /3 and 1/6 filling, 
respectively (and their respective particle-hole values). 

In the following we will be interested in the strongly 
anisotropic limit Jz \Jxy\ where we will see that the 
dimer physics emerges. Since the diagonal part in (1) can 
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FIG. 2. (Color online) Ground state magnetization curves 
for a system of size L = 6 and different values of J^y, in the 
ground state {Jz = 1 is taken as the unit of energy). In¬ 
set: Magnetization curve for Jz = 0 (Ising limit). The three 
plateaus at magnetizations m = 0,1/6,1/3 are clearly visi¬ 
ble for the lowest value of the off-diagonal coupling constant 
Jxy = —0.06, and gradually disappear as gets larger. 

be written as: 

= ( 2 ) 

O O 

with Nq the number of hexagons, every state satisfying 
the constraint Sq = 0,1, 2,3 (depending on the value of 
h) on all the plaquettes is a ground state of the classi¬ 
cal part Hq. First excited states are separated by a gap 
of magnitude Jz- As a consequence, the magnetization 
curve m{h) (with normalization m = /N) displays 

plateaus of width Jz at the values m = 0,1/6,1/3, (see 
the inset of Fig. 2) and saturation at 1/2 (m = (3 —n^)/6 
where = 3, 2,1,0 is the number of down spins per 
hexagon) which are expected to survive for a finite value 
of Jxy As an illustration, we present in Fig. 2 the magne¬ 
tization curve in the ground state for a moderate system 
size L = 6 and finite small values of Jxy, as obtained from 
QMC simulations (see Sec. IV for details and a finite-size 
scaling analysis of the size of the m = 1/6 plateau). 

The above construction leads to magnetization 
plateaus where the ground states degeneracies scale ex¬ 
ponentially with the system size [10]. Dynamics induced 
by the off-diagonal terms will lift this degeneracy, as pre¬ 
sented below. 


B. Mapping to generalized quantum dimer models 

A simple way to visualize the local constraints S'q = 
0,1, 2,3 (for all hexagons) in the Ising limit is to faithfully 
map the ground state configurations to those of dimers: 


FIG. 3. (Golor online) Representation of one of the ground 
states of model (1) for J^y = 0 and m = 1/6 in terms of 
dimers (red) on the underlying triangular lattice (blue). Loop 
configurations are formed due to the constraint of two dimers 
per triangular lattice site. 

draw a dimer on a link of the triangular lattice formed 
by the centers of hexagons when the middle of the link is 
occupied by a down spin, and leave the same link empty 
for an up spin. The magnetization plateaus m = (3 — 
n^)/6 are then reproduced by having a constraint such 
that every triangular lattice site must be touched by 
dimers. This procedure is illustrated in Fig. 3 for the 
m = 1/6 (n^ = 2) plateau with two dimers per site, which 
allows us to visualize the ground state spin configurations 
of the diagonal part Hq as fully-packed self-avoiding loops 
on the triangular lattice. 

The generalized quantum dimer models appear when 
treating Hxy in degenerate perturbation theory. Clearly, 
first-order processes do not contribute, since applying 
the spin-flip term only once produces configurations with 
two excited hexagons, thus outside of the ground state 
manifold. In second-order, two types of processes are 

allowed. The first process simply flips the same pair of 

f2l 

spins twice, merely shifting the total energy by = 

—9{l—m‘^)NQj^y/Jz. The second process involves a flip- 
pable bow-tie plaquette, as pictured in green in Fig. 4. 
The two pairs of antiparallel spins are flipped, generating 
another state fulfilling the ground state constraint. The 
second-order effective Hamiltonian is thus given by the 
ring-exchange Hamiltonian [10]: 

= -Jring * 53 +^ 4 - + h.c.), (3) 

where the sum runs over all flippable bow-tie plaquettes 
(labels 1, 2, 3 and 4 denote the exterior sites of the bow- 
ties) and the ring-exchange amplitude is given by Jring = 
^Jxy! Jz- Quite importantly, the sign of Jxy is immaterial 
which allows numerical simulations of the original micro¬ 
scopic Hamiltonian (1) with QMC techniques, by consid¬ 
ering the model with a ferromagnetic Jxy < 0 coupling. 
Also in order to obtain the ring-exchange Hamiltonian 
(3), only the XY-interactions on the nearest-neighbor 
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FIG. 4. (Color online) Off-diagonal second-order process of 
Hxy on a flippable bow-tie plaquette (green) for n 4 , = 3. The 
antiparallel spins of two flippable pairs are exchanged, corre¬ 
sponding to a flip of two parallel dimers on the underlying 
triangular lattice (blue). 


{J^y) and the next-nearest-neighbor (J^'^^) bonds are 
crucial: 


'.^rinp' 


j. 


( 4 ) 


From this, one can see, as stated previously, that sup¬ 
pressing the further-neighbor XY interactions in H^y 
[Eq. (1)] does not alter physics qualitatively and simply 
rescales Jring by a factor 1/2. We exploit this property 
in Sec. IV. 

In the dimer language, the ring move on a bow-tie of 
the original lattice corresponds to the flip of two paral¬ 
lel dimers on a diamond of the dual triangular lattice 
(Fig. 4). This is nothing but the off-diagonal term of 
the quantum dimer model [7] on a triangular lattice [8], 
which, as usual, is conveniently extended with a diagonal 
potential (V) term which counts the number of flippable 
plaquettes (diamonds): 




i:( 

A7 


A7 


/)(2:i 


( 5 ) 


In the Ising limit Jz ^ Jxy, low-energy physics of the 
BFG model is expected to be described by this Hamil¬ 
tonian for V = 0. The following four-spin interaction 
[10, 13] 


Q45] {(1/2 - St ) (1/2 + St ) (1/2 - 5 |) (1/2 + St ) 

[X 

+ ( 1/2 + ^/) ( 1/2 - St) ( 1/2 + 5 |) ( 1/2 - 5 |)} 

( 6 ) 


(1,2 ,3,4 label the four spins on a bow-tie in a anti¬ 
clockwise way) reduces, in the same limit, to the V-term 
(V — Q) as it counts the number of flippable plaquettes. 
For the m = 1/3 plateau corresponding to one dimer 
(or, one { spin) per hexagon, one may use instead the 
inter-hexagon (two-spin) third-neighbor interaction 


E 


qz qz 
Dj^ Oj 


inter-hexagon 


( 7 ) 


that is not contained in model (1) to realize the dimer 
interaction V = J^\ 

Thus, one sees that the effective model is identical for 
the three magnetization plateaux, the only difference com¬ 
ing from the constraint on the Hilbert space. For the con¬ 
straint of one dimer per site (corresponding to the plateau 
m = 1/3), we have the usual QDM [8], while in the case 
with two dimers per site (the plateau at m = 1/6), a 
QLM is obtained as announced in the introduction. The 
original model derived by BFG [10] describing the m = 0 
plateau corresponds to three dimers per site. Table I 
summarizes the nature of the effective models for the 
different plateaus. 

The mapping to dimer models [including the diagonal 
V-term in Eq. (5)] allows us to benefit from the accumu¬ 
lated knowledge on this family of models [7-9, 23-30]. 
QDMs generically admit a Rokhsar-Kivelson [7] (RK) 
point V = t, where the ground-state wave-function is 
exactly known. For the QDM on the triangular lattice, 
the ground state is exactly shown to be a topological 
gapped Z 2 liquid phase [8, 30, 31], which furthermore 
extends in the region V/t < 1. In the three-dimer model, 
BFG confirmed the presence of a topological phase at 
the RK point [10], showing in particular that visons are 
gapped, implying the presence of deconfined fractional¬ 
ized spinons [32]. The authors of Ref. 10 also specu¬ 
lated that the Z 2 liquid phase could extend beyond the 
V = 0 point corresponding to the original microscopic 
spin model, in contrast to the one-dimer case [8, 25] for 
which the liquid phase ends at {V/t)c — 0.8 [25]. An 
exact diagonalization study [13] of the Hamiltonian (3) 
shows the presence of a vison gap in a regime including 
the BFG point (V = 0), thereby supporting the above 
suggestion. However, the presence of the spin liquid was 
only firmly established using QMG simulations on large 
systems. First, a complete numerical phase diagram of 
the model (1) with ferromagnetic Jxy and h = 0 was ob¬ 
tained [14]: at strongly negative Jxy!Jz, a planar (ie., 
superfluid in the bosonic language) phase accompanied 
by the breaking of the U(l) symmetry is found, as ex¬ 
pected. When the magnitude of the spin-flip term was 
increased, a continuous phase transition to an apparently 
featureless insulating phase was observed. The phase di¬ 
agram was also extended to finite temperature [33]. Fi¬ 
nally, a direct evidence of the topological nature of the 
insulator is obtained by computing the topological en¬ 
tanglement (Renyi) entropy [16], which was evaluated 
as log 2 as expected for a Z 2 spin liquid phase. As a 
consequence of the condensation of fractionalized excita¬ 
tions, the transition between the planar phase and the 
featureless insulator is an exotic 3D XY* quantum criti¬ 
cal point [15, 34]. 

While the nature of the phases encountered at the m = 
0 and m = 1/3 plateaus is well understood from this 
previous set of study, the ground state physics on the 
m = 1/6 plateau has never been investigated to our best 
knowledge. In the following, we focus on this situation 
by performing numerical simulations first on the effective 
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FIG. 5. (Color online) Phase diagram of the QLM Eq. (5) as a function of Vjt, containing the three different phases discussed 
in Sec. Ill A. The boundary between stripe phase and Z 2 spin liquid is determined using topological degeneracy and stripe 
structure factor. 


TABLE I. Three plateaus of S' = 1/2 BEG model and the 
corresponding effective models. 


plateaus q{S — m) 

number of dimers/site 

effective model 

m — 0 3/2 

three 

generalized QDM 

m = 1/6 1 

two 

QLM 

m = l/3 1/2 

one 

QDM 


QLM (Sec. Ill), and then on the microscopic spin model 
(Sec. IV). 


III. PHASE DIAGRAM OF THE QUANTUM 
LOOP MODEL 

We consider here the effective model (5) in the case 
of two dimers per site, where allowed states are repre¬ 
sented by configurations of self-avoiding loops such as 
represented in Fig. 3. Motivated by the original spin 
model Eq. (1) at the m = 1/6 plateau, we are primarily 
interested in the nature of the ground state of the effec¬ 
tive QLM at V/t = 0, however we will also investigate the 
nature of the surrounding phases in the phase diagram. 
We first give, in Sec. Ill A, general arguments on the 
structure of the phase diagram as well as on the topolog¬ 
ical properties of the loop configurations, before comple¬ 
menting this analysis in Sec. Ill B with QMC simulations 
of the loop model. Anticipating the results obtained, we 
present in Fig. 5 the ground-state phase diagram of the 
QLM on the triangular lattice to illustrate the following 
discussion. 


A. General considerations of the phase diagram of 
the quantum loop model 

When Vjt —>■ — 00 , the ground-state energy is mini¬ 
mized by configurations where loops form straight lines 
along one of the three lattice directions, as shown in the 
left part of the phase diagram presented in Fig. 5. As 



FIG. 6 . (Golor online) Top left: Triangular clusters used in 
the QLM study. Bottom left: notations for the loop segment 
occupation number. Right: First Brillouin zone of the tri¬ 
angular lattice, with the reciprocal space vectors hi and 62 . 
The high symmetry points K = ( 47 r/ 3 , 0), M = (vr, 7r/-\/3) 
and A = (tt, 0) required for the non-flippable states are rep¬ 
resented. 


aligned up and down spins form alternating stripes in the 
corresponding configurations of the original spin mod¬ 
els (see the right panel of Fig. 12), we call this a stripe 
phase. In this gapped ordered phase, some of the rota¬ 
tions and reflections of the triangular lattice are sponta¬ 
neously broken {i.e., the point group changes from Cqv to 
C' 2 v), leading to the three-fold degenerate ground states 
(in the QDM, the corresponding columnar phase further 
breaks translation symmetry). 

At the RK point V/t = 1, the ground state is given by 
the equal-weighted sum of all fully-packed loop configu¬ 
rations on the triangular lattice. To our best knowledge, 
the equivalent classical problem has not been studied be¬ 
fore. As the triangular lattice is non-bipartite, we expect 
that the loop segment correlations are short-ranged, in¬ 
dicating a liquid phase presumably with a gap, as also 
found for the RK point of the QDM [8] . 

When V/t > I, the ground states are readily found 
to have zero energy and correspond to non-flippable loop 
configurations {i.e. for which the kinetic term vanishes), 
which is again similar to the staggered phase of the QDM. 
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While the construction of those states is quite straightfor¬ 
ward in the QDM, we have found that for L x L systems, 
the QLM accepts two different families of non-flippable 
configurations, both of which are shown in the right part 
of the phase diagram of Fig. 5. The first family is ob¬ 
tained by starting from the staggered states of the QDM 
[8, 23], which are 12-fold degenerate. On the left, one of 
those 12 states is represented in red, with only two of the 
three types of links (corresponding to the three lattice di¬ 
rections 0, 1, and 2 in Fig. 6) occupied. A non-flippable 
loop configuration can then be generated by adding par¬ 
allel (blue) dimers on half of the remaining type of links, 
for which there are two possibilities (a one-step transla¬ 
tion in the direction of these parallel dimers also creates a 
different non-flippable state). The same prescription can 
be applied starting from the other QDM staggered states, 
leading to a total degeneracy of 2 x 12 = 24. The second 
family is obtained by tiling the lattice with the shortest 
possible triangular loops, as depicted on the right part of 
Fig. 5. The degeneracy of this family is readily found to 
be 6. Considering the Brillouin zone drawn in the right 
panel of Fig. 6, the first family of non-flippable states 
requires the cluster to have both A and M points (i.e. L 
multiple of 4), while the second one needs the K point 
(i.e. L multiple of 3). We finally remark that the na¬ 
ture of non-flippable states (or absence thereof) can be 
different for other types of clusters (rectangular clusters 
for instance). 

From these considerations, the point V = Q lies be¬ 
tween the two limiting cases Vjt ^ —oo (stripe) and 
Vft —)■ 1~ (liquid). As both the stripe and liquid phases 
are gapped, they should extend in a finite region of phase 
space close to these limits, and it is not clear a priori in 
which phase the point D = 0 is located. Of course, it 
is also possible that one or several intervening crystalline 
phases arise: natural candidates are equivalent of plaque- 
tte, or phases observed in QDMs [8, 23, 25, 35]. 

This can only be determined by exact numerical calcula¬ 
tions, which will be presented below in Sec. IIIB. 

It is interesting to note at this stage that on a mani¬ 
fold with non-zero genus (cylinder, torus, etc), the QLM 
model possesses topological sectors which can be defined 
in the same way as in the QDM case. As illustrated 
in Fig. 7 for a torus, the parity of the number of loop 
segments crossing a dashed line in either the x or y di¬ 
rections of the torus is conserved by the dynamics of the 
Hamiltonian. On the torus, this defines a pair of Z 2 in¬ 
variants Pa; = ±1 and Py = ±1, and thus divides the 
full Hilbert space into four distinct topological sectors, 
which we label by {px,Py)- The existence of these four 
topological sectors is the key to identify the Z 2 spin liquid 
phase close to V/t = 1, as e.g. in the QDM on the tri¬ 
angular lattice [8] . One may wonder if the Z 2 spin liquid 
found here is the same as that in the usual triangular- 
lattice QDM since even (odd) number of dimers exiting 
from each site in the QLM (QDM) suggests a mapping 
to an even (odd) Ising gauge theory [36]. However, this 
difference is not important in considering the underlying 


Pa: = +1 / 



FIG. 7. (Color online) Drawing two lines on the triangular 
lattice allows to define two conserved Z 2 invariants {px,Py) as¬ 
sociated to the parity of the number of loop segments crossed 
by the line. In this example, the number of loop segments 
pointing in the vertical y (horizontal x) direction which cross 
the line is odd (even), corresponding to Py = —1 (px — 1). 
This number does not depend on the precise location of the 
line (parallel lines shifted by one or several lattice units result 
in the same parity). The flipping terms of the Hamiltonian, 
acting for instance on the shaded plaquette (other flips are 
possible in this configuration), do not change these two pari¬ 
ties. Note that a loop segment located on a diagonal bond of 
a plaquette contributes to both px and py. 


nature of the spin liquids. In fact, by explicitly construct¬ 
ing the ground state of the toric code [37] on a triangular 
lattice [38], one can show that both spin liquids share 
the same quasiparticle (anyon) contents and exhibit the 
same ground-state degeneracy originating from the Z 2 
flux, which we will use in the following. 

In the topological Z 2 spin liquid phase, the four topo¬ 
logical sectors (on a torus) corresponds to the four differ¬ 
ent ways of threading the Z 2 flux through the two periods 
of the system that cannot be eliminated by finite-order 
perturbations [39]. Therefore, the four sectors, in the 
topological spin liquid phase, are separated from each 
other by exponentially small gaps {topological gaps) and 
get degenerate in the thermodynamic limit. As was pre¬ 
viously observed for the QDM [23], we expect the liquid 
phase to be detected by the closing of the above topolog¬ 
ical gaps. 

On the other hand, in the ordered stripe phase, topo¬ 
logical sectors which do not contain the pure stripe con¬ 
figurations (shown in the left part of Fig. 5) cannot op¬ 
timize the attractive D-term and should have a higher 
energy: correspondingly, the topological gap should be 
finite even in the thermodynamic limit, providing us with 
a practical method to distinguish the two phases. One 
should however be careful in identifying the sectors to 
which the stripe configurations belong. Indeed on a torus 
with Lx = Ly = L, the QLM can be studied for even and 
odd L (in contrast to the QDM that can only host even- 
L samples). When L is odd, the three stripe configura¬ 
tions shown in Fig. 5 are located in the three different 
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V/t V/t 

FIG. 8. (Color online) Topological gap At as a function of 
Vjt, for different sample sizes L. The finite-size effects are 
different for even left panel) and odd (At‘^‘^; right 

panel) values of L. In the thermodynamic limit, we expect 
that in both cases, At if finite for the stripe phase Vjt < 
{yit)c, and vanishes for the Z 2 spin liquid phase (F/t)c < 
V < 1. Left inset: Finite-size scaling as a function of l/L 
for the even-L samples provides an estimate {V/t)c — 0.70(5) 
of the value where the gap first vanishes. Right inset: Zoom 
close to the RK point for odd-L samples. 


sectors (—1,1), (1,-1) and (-1,-1) which are strictly 
identical (the number of configurations is the same and 
the corresponding blocks of the Hamiltonian are iden¬ 
tical): we define the topological gap for odd-L samples 
as = Eo{l,l) - i;o(-l,-l), where Eo{px,Py) is 

the ground state energy in the {px,Py) sector. For even 
L, the three stripe configurations belong to the (1,1) 
sector, and all other sectors are identical, and we use 

a-®" = f;o(-i,-i)-f;o(i,i). 


B. Phase diagram: QMC results 

We now turn to QMC simulations of the QLM on 
a torus (for clusters with the geometry represented in 
Fig. 6, with Lx = Ly = L and N = L'^ sites), as de¬ 
fined in Eq. (5). We use the reptation QMC [40], which 
projects a given initial configuration belonging to each 
topological sector onto the ground state in the same sec¬ 
tor. We use a formulation of reptation QMC similar to 
the one presented in Ref. 35: to accelerate the conver¬ 
gence, we work in a continuous-time formalism as well as 
with a guiding wave function that provides a finite fugac- 
ity for each flippable plaquette of the triangular lattice 
(an optimized value for the fugacity is obtained from a 
short preliminary run). We simulated even-L samples up 
to L = 10 (and one data point at L = 12), and odd-L 
samples up to L = 9. 

We first present the results for the topological gap At, 
as obtained by computing the ground state energy in each 
non-equivalent sector using the standard mixed estima¬ 


tor [40]. Fig. 8 shows At as a function of V/t for even 
(left panel) and odd (right panel) L, for different system 
sizes. For even samples, the topological gap A^™" has 
a monotonous behavior: it is clearly finite for low val¬ 
ues of V/t and vanishes with system size close enough 
to the RK point V/t = 1. Finite-size scaling as a func¬ 
tion of 1/L (left inset of Fig. 8) indicates that the gap 
vanishes around {V/t)c — 0.70(5), separating the stripe 
phase for V/t < {V/t)c from the topological Z 2 liquid 
phase for {V/t)c < V/t < 1. For the largest samples, our 
data at V/t = 0.8 and V/t = 0.9 show that the ground 
state energies in different topological sectors are identi¬ 
cal within error bars (we cannot resolve the exponentially 
small splitting in the liquid phase). 

The gap for the odd-L samples displays an appar¬ 
ently different behavior, with the topological sector (1,1) 
having a higher energy for low enough K/t, resulting in 
At'^'^ > 0 as the K-term favors the sectors (—1,1), (1,-1) 
and (-1,-1). On the other hand, when the system is 
close enough to the RK point, the sector (1,1) has lower 
energy (resulting in AJjd'^ < 0 in our definition of ASjd'^) 
for finite-size samples. If we assume that the ordering of 
the four levels in the spin liquid phase is the same as that 
for the even-L samples, the sign change marks the tran¬ 
sition. The value of V/t at which = 0 varies with 
system size and appears to extrapolate to a consistent 
value for the critical point {V/t)c — 0.70(5) (see right 
inset of Fig. 8). Here again, for the largest odd-L sam¬ 
ples, the energies in all topological sectors are indistin¬ 
guishable within error bars for V/t = 0.8 and V/t = 0.9, 
indicating the topological degeneracy. 

The topological gap for both odd and even sample is 
clearly finite for V/t = 0 (with a value ~ 1.3t, as esti¬ 
mated from the even-L data, for which there is almost 
no size dependence), which indicates that this point is 
located in an ordered phase (presumably, in the stripe 
phase). In order to confirm this, we computed the loop 
segment occupations (using middle-slices in the reptile 
representation [40]) in the topological sector (1,1) (which 
hosts all three stripe configurations) for the even-L sam¬ 
ples. This allows to determine the stripe structure factor, 
which we define as: 

-S'stripe = ^ ^ 

a,/3=0,1.2 

= ^ E (iv^)-E(^«^/5) • 

ya=0,1.2 J 

where ri^Q, = 1 if there is a loop segment at site i occu¬ 
pying a link in the direction a{= 0,1,2; the three values 
taken by a correspond to the three different directions as 
represented in the left bottom part of Fig. 6) and rii^a = 0 
otherwise. In the second line of Eq. (8), we have intro¬ 
duced the notation Na = for the total number of 

bonds in the direction a. The structure factor divided by 
the sample size is expected to converge to the square of 
the stripe order parameter rustripe in the thermodynamic 




















limit: ^stripeZ-ZV with a non-zero rngtripe if the 

ground state exhibits the long-range stripe order. 



FIG. 9. (Color online) The stripe structure factor S'stripe as a 
function of inverse system size, for different values oiVjt. 

The stripe structure factor S'stripe/fV is shown in Fig. 9 
for different values of the coupling Vjt as &. function of in¬ 
verse size 1/N (for L = 4,6,8). We expect that Sstn-pe/N 
should saturate in an exponential way (due to the finite 
correlation length) to > 0 in the stripe phase, 

and to 0 in a phase with no stripe order (at the criti¬ 
cal point, it should decay to 0 with a non-trivial power). 
Due to the moderate system sizes accessible, we cannot 
expect to see the exponential decay: the extrapolation 
to the thermodynamic limit is indeed delicate for large 
values of V/t. Considering the data for V/t < 0.6 (see 
Fig. 9), we observe that the structure factor extrapolates 
to a non-zero value. From this data set, we obtained an 
estimate for the critical point {V/t)c = 0.65(15), which is 
less precise, though consistent, than the one determined 
above by the topological gap. 

Note that, in analogy with the situation in the QDM, 
one may suspect that other types of symmetry-breaking 
order (such as plaquette order) could also exhibit a non¬ 
zero value of TTistripe- To clarify this, it is useful to con¬ 
sider the histograms of the occupations of loop segments 
appearing in the Monte Carlo simulations, as was done in 
similar simulations for valence bond crystals [41-45]. To 
have a two-dimensional representation of the histograms 
for the loop segments occupying one of the three lattice 
directions, we introduce = ^(2iVo — Ni — N 2 ) and 

Hy = 2 ^(—fVi -I- N 2 ) (with Na being the total number 
of bonds in the direction a). The three line configura¬ 
tions shown in the left part of Fig. 5 (having respectively 
Nq, Ni or N 2 equal to N) occupy, in the (iJa,, Hy) repre¬ 
sentation, the three vertices (1, 0), (—1/2, ±-\/3/2) of an 
equilateral triangle in which the histogram is enclosed, as 
represented in Fig. 10. The left inset of Fig. 10 presents 
the histogram obtained for a moderate system size L = 6 
and V/t = —0.3, which is already located deep inside 



FIG. 10. (Golor online) Expectation value W 3 of the quanti¬ 
fier of Z 3 anisotropy of histograms of loop occupations, as a 
function oiV/t for three even lattices sizes L = 4, 6, 8. Insets: 
Golor-coded normalized histograms of loop occupations num¬ 
bers (using a color interpolation scheme), deep in the stripe 
phase (left inset, L = 6, V/t = —0.3) and closer to the crit¬ 
ical point (right inset, L = 6, V/t = 0.3). Vertical axis of 
histograms is Hy (ranging from —%/3/2 to \/3/2), horizontal 
axis is Hx (from —1/2 to 1). The histograms are enclosed in 
an equilateral triangle pictured close to the right inset. 


the stripe phase: we clearly observe a Zs-symmetric fea¬ 
ture with three peaks (corresponding to high occupa¬ 
tions) close to the corners of the triangle corresponding 
to the three line configurations. The distance of these 
peaks from the origin (0, 0) gives a measure of the order 
parameter Wstripe- 

On the other hand, the angles 0 = avctan^Hy/H^) as¬ 
sociated to the three peaks in this two-dimensional rep¬ 
resentation also carry important information on the pre¬ 
cise nature of the phase. For instance, the angles associ¬ 
ated with the three stripe patterns are 0, 27r/3 and 47r/3. 
Other candidate phases (such as the equivalent of pla¬ 
quette phases) would show up in the histograms as sharp 
peaks at different angles. We find that, as V/t gets closer 
to the critical point {V/t)c, the typical extents (radius) 
of the histograms get smaller (as expected from the de¬ 
crease of the order parameter), and that at the same time 
the shape of the histograms change from a Za-symmetric 
one to a circular-symmetric one as is exemplified in the 
right inset of Fig. 10. This is reminiscent of the C/(l) 
valence bond histograms observed in the vicinity of de- 
confined quantum critical points [41, 42, 44, 46] or even 
at moderate proximity of the RK point of the QDM on 
the square lattice [47, 48]. 

To quantify this effect, we calculated a measure of the 
Za-symmetry of the histograms Wa = (cos(30)) using for¬ 
mulations similar to those in Refs. [42, 44-46, 49]. For a 
pure stripe phase, we have Wa = 1, while for a circular 
histogram W 3 = 0. The variation of W 3 as a function V/t 
with different system sizes displayed in the main panel of 
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Fig. 10 confirms the behavior observed visually in the 
histograms: for the region V/t < 0.1, W 3 tends to ap¬ 
proach unity with increasing system sizes, while in the 
region 0.1 < V/t < 1, W 3 appears to vanish within the 
system sizes available to us. While therefore we cannot 
exclude an intermediate different crystalline phase, it ap¬ 
pears unlikely as there is no direct evidence either in the 
energy or in the structure factor of a phase transition (see 
a similar situation for the square lattice QDM [47, 48]). 
Thus, we expect that on larger systems W 3 would also 
tend to approach unity for V/t < {V/t)^. Note that the 
computations of W 3 are difficult as we need long Monte 
Carlo runs to be fully ergodic (i.e., in order not to get 
locked in one type of line configurations, in particular in 
the stripe phase), and also to ensure that statistical fluc¬ 
tuations are not too important when the order parameter 
is small. We may speculate that these circular histograms 
are indications of the Z 3 anisotropy being a dangerously 
irrelevant variable at the critical point {V/t)^. Although 
the precise characterization of the transition between the 
Z 2 spin liquid and the stripe phase is an interesting issue, 
it is beyond the scope of the present work. At the point 
V/t = 0 oi our main interest, our QMC data (see Figs. 9 
and 10 ) clearly indicate that the ground state of the ef¬ 
fective Hamiltonian (3) is located in the stripe phase. 
The resulting phase diagram obtained from these QMC 
results as well as the general considerations in Sec. Ill A 
is given in Fig 5. 


IV. NUMERICAL SIMULATIONS OF THE 
MICROSCOPIC MODEL 

In this section, we present the results of the QMC sim¬ 
ulations for the original spin-1/2 Hamiltonian (1) at the 
m = 1/6 plateau. We used the Stochastic Series Expan¬ 
sion (SSE) algorithm [50, 51] combined with a plaquette 
generalization [52] that is necessary to circumvent the 
freezing encountered in the bond formulation of the SSE 
when Jz is large [53, 54]. In addition to this difficulty, 
capturing ground state physics requires to reach a very 
low temperature-range characterized by the small energy 
scale T ~ J^ing ^ Jxy/Jz of the effective model. Never¬ 
theless, we have been able to simulate LxL systems, with 
N = 3L^ sites, up to L ~ 12, or slightly beyond L ~ 18 
for some observables. Eor numerical ease, we restrict 
ourselves to the truncated model where only the nearest- 
neighbor interactions (i.e., the black bonds in Fig. 1) are 
retained in Hxy (see the discussions in Sec. HB). 

In the following, we set = 1 s.s the energy scale, 
and use ferromagnetic values Jxy < 0 (for the QMC 
simulations to have no sign problem) and h = 1.03 for 
the magnetic field (corresponding to the middle of the 
TO = 1/6 plateau; see Fig. 2). In order to locate the 
transition point out of the plateau phase as a function 
of Jxy, we performed a finite-size analysis of its width. 
Resulting data are shown in Fig. 11, and the transition 
point is estimated to be around Jxy — —0.081. For val- 



FIG. 11. (Color online) Finite-size scaling of the m = 1/6 
plateau width. The plateau is found to disappear around 

Jxy — —0.081. 


ues Jxy < —0.081, the system is thus in a planar phase 
(superfluid in the bosonic language), and in the follow¬ 
ing we focus on the plateau region —0.081 < Jxy < 0 . 
In Sec IV A, we examine the nature of the plateau phase, 
and discuss what could be the scenarios for the transition 
to the planar phase. Sec IV B is devoted to an alternative 
characterization of the phase using the Renyi entangle¬ 
ment entropy. 


A. Nature of the m = 1/6 plateau 

For magnetization plateaus, as has been mentioned in 
the introduction, it is useful to consider the quantity 
q{S — to), where q is the number of sites per unit cell 
(g = 3 for the Kagome lattice). Indeed, it has been 
shown analytically [17, 19] that when this quantity is 
fractional, the ground state on this plateau must be de¬ 
generate, i.e., either in a crystalline state with magnetic 
superstructures or in a topological state, while a non¬ 
degenerate featureless state with a gap is forbidden. If, 
on the other hand, it takes an integer value, which is the 
case here at the to = 1/6 plateau (see Table I), all these 
possibilities are allowed. 

Guided by the results of the effective QLM in Sec HI, 
we expect the plateau phase to be in a ordered phase 
where the down spins form stripes. This is indeed evident 
when plotting the connected S^-S^ spin correlation func¬ 
tion in real space, as shown in the left panel of Fig. 12, 
where one readily sees the presence of stripes. The cor¬ 
responding phase breaks the six-fold rotation symmetry, 
and is three-fold degenerate. Cartoon pictures of the spin 
and dimer configurations in this phase are drawn on the 
right panel of Fig. 12. 

In order to confirm the long-range nature of this order, 
we compute the Fourier transform of the above spin-spin 
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FIG. 12. (Color online) Left: Connected correlation function(S'f)c = {Si Sj) — {Si){Sj) on the m = Ijd plateau (for 
Jxy = —0.078, T = Jring/8 On a 1/ = 8 cluster with periodic boundary conditions). Reference site is indicated as a black 
circle. Blue (respectively red) denote positive (resp. negative) values, and the radius is proportional to the absolute value of 
the correlation. Data close to the reference site have not been plotted for readability. Right: Cartoon representation of the 
spin and dimer configurations in one state of the stripe phase, with the Kagome lattice in black and the underlying triangular 
lattice in blue. 


correlations, i.e., the diagonal spin structure factor: 

3,k 

where is the position of spin N = 3L^, q = {qx,qy) 
and the average local magnetization is (Sj) = m = 
1/6. For an ordered phase at wave-vector Q, S'(Q) di¬ 
verges as N, and for the three-fold degenerate solid, the 
three Bragg peaks are located at wave-vectors Qi = 
(0,47r/-\/3), Q 2 = (27r, 27r/-\/3) and Q 3 = {2tt,—2-k/ y/Z). 
Note that this solid can be accommodated on any cluster 
used for the simulations. Finite-size scaling of S'(Q)/Ai as 
a function of 1 /L for several values of Jxy in the plateau 
phase clearly indicates the presence of the stripe phase 
(see Fig. 13). In agreement with the QLM study at 
V/t = 0, the order develops already for small system 
sizes. 

Two comments are in order here. First, 5'(Q) used in 
Fig. 13 is the average of the values at the three wave- 
vectors Qq, which is valid when the simulation is per¬ 
fectly ergodic. However, as mentioned earlier, the large 
anisotropy Jz j Jxy 13 in the parameter range of interest 
can create difficulties in the QMC simulations [54]. By 
monitoring the three values iS'(Qa) separately, we have 
checked that the simulations are ergodic for the sizes 
L < 9 considered here. However, for smaller values of 
\Jxy\lJz or larger L, we indeed noticed a quick loss of 
ergodicity, even with the plaquette algorithm, translat¬ 
ing into very different values at the three wave-vectors. 
The same behavior was observed in Ref. 53 which used a 
similar QMC algorithm. Second, we sometimes found a 
non-monotonic behavior of the structure factor S{Ci)/N 
for larger sizes (see Fig. 13). Again, a similar behav¬ 
ior was reported in Ref. 53 where the increase of the 
order parameter was interpreted as the threshold where 
stripes start developing in the system. In practice, this 



FIG. 13. (Color online) Finite-size scaling of the (scaled) 
structure factor S{Q)/N inside the magnetization plateau. 
Despite the relatively small available sizes, we can reach the 
regime of saturation, confirming the striped solid anticipated 
from the QLM. 


complicates a proper extrapolation of the order param¬ 
eter to its thermodynamic limit value. We can how¬ 
ever compare the order of magnitude obtained with the 
value for a perfectly equal superposition of the three 
ordered states, for which a straightforward calculation 
yields S{Q)/N = 8/81 ~ 0.0987. 

We conclude this section by a discussion on the possi¬ 
ble scenarios for the quantum phase transition between 
the plateau (stripe) crystal and planar phases. Since the 
stripe crystal breaks discrete lattice symmetries, one ex¬ 
pects, on general grounds, that the transition will be of 
first order type, or that both phases coexist (in a small 
region of the phase diagram). This latter scenario implies 
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the existence of an intervening supersolid phase, whose 
presence was found in several similar models on other 
frustrated lattices [55-57]. Due to the limited systems 
sizes available to us, we have not been able to favor one 
of these two scenarios. We observed (data not shown) an 
apparent crossing in the spin stiffness, with critical ex¬ 
ponents compatible with a continuous phase transition, 
and no sign of a double peak structure in the kinetic en¬ 
ergy histograms. However, it has been shown in similar 
situations that the correct nature of the phase transition 
may be difficult to capture. In the {Jxy,h) plane, varying 
the XY interaction J^y at the constant field h = 1.03 
corresponds to a transition point located close to the tip 
of plateau phase (i.e. the tip of the insulator lobe in the 
bosonic picture), for which previous studies have revealed 
that such a crystal-planar phase transition may appear 
to be continuous while it is in fact weakly first-order [58- 
60]. We believe that such a picture should also apply 
here, although proving this definitively would imply sim¬ 
ulations on very large samples, which is currently out of 
reach for the present model. 


B. Renyi entropy in the crystalline phase 

We have shown, using a conventional approach (be. by 
guessing a broken symmetry and then measuring the cor¬ 
responding structure factor), that the m = 1/6 plateau 
corresponds to a stripe crystal shown in Fig. 12. It is in¬ 
teresting to check whether other means could detect the 
nature of the ground state without a priori knowledge 
on the broken symmetry. An elegant approach to obtain 
the dominant order parameter in an unbiased way is pro¬ 
vided by the correlation density matrix [61], but it is not 
easily accessible within QMC simulations. Instead, we 
will focus here on using the entanglement entropy as a 
means to access the nature of the ground state. 

Indeed, in the recent years, a large number of studies 
have shown how the scaling behavior of a block entangle¬ 
ment entropy can give access to some ground state prop¬ 
erties, such as the central charge in one-dimension [62] 
or the topological order underlying a given ground state 
wave function [63, 64]. Generically, any entanglement en¬ 
tropy (whether Renyi or von Neumann) will scale with 
the “area” of the boundary between the block and the 
rest of the system (the so-called “area-law”; see Ref. 65 
for a review), possibly with interesting (universal) sub¬ 
leading terms: 

Sq{A) = aq(.A+ dq ( 10 ) 

where q is the Renyi index, £a the length of the bound¬ 
ary of block A and dq a constant for instance. Topo¬ 
logical phases can be detected using the constructions of 
Levin-Wen [63] or Kitaev-Preskill [64] (see also Ref. 66 ), 
that enables one to extract a negative constant term 
dq = —7 < 0, independent of q [67], where 7 is the topo¬ 
logical entanglement entropy. This was used for instance 
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FIG. 14. (Color online) Renyi block entanglement entropies 
Sq{A) {q — 2, 3) on the m = IjG plateau, as a function of L, on 
various systems with N = 3L^ sites. Simulations parameters 
are J^y = —0.078, T = Jring/4 and h = 1.03. In each case, 
the A block corresponds to a half-system cylinder containing 
(L/2) X L unit cells (see inset). The positive correction log(3) 
is shown by a star. 

to conclude that the ground state of the BFG model at 
m = 0 is indeed a Z 2 spin liquid with 7 = log 2 [16]. 

The constructions of Refs. [63, 64] use a subtraction of 
terms which all scale like the area of the subsystems con¬ 
sidered, which can be troublesome in QMG simulations 
which inevitably exhibit finite error bars. Gomputing 
the scaling of Sq{A) for a single block A (without using 
subtraction schemes), is in general computationally sim¬ 
pler. This is what we performed for the case of our dis¬ 
crete symmetry breaking (Cev '—t C' 2 v) on the m = 1/6 
plateau, taking for A the geometry of a cylinder block 
embedded in a torus sample (see the inset of Fig. 14). 

We expect that for such a symmetry-broken state, the 
sub-leading term will be a positive constant dq = logg, 
where g is the number of degenerate ground states (note 
the sign of dq; see Appendix A). We first remark that, 
for the simple type of order on the m = lj6 plateau, 
the precise block geometry is not very important, which 
allows us to consider the convenient cylinder block ge¬ 
ometry. Second, it is quite important for this result to 
hold that all g degenerate states (quasi degenerate on a 
finite system but with an exponentially small splitting) 
are observed in the QMG simulations. Hence, ergodicity 
must be ensured to extract the degeneracy, which we ob¬ 
serve is the case if we are not too deep inside the stripe 
phase and system sizes are not too large. Note that such 
a positive correction gives a vanishing contribution to the 
previously mentioned geometric constructions [63, 64], as 
they are precisely built to solely capture long-range en¬ 
tanglement. Our numerical data are presented in Fig. 14. 
We have computed the Renyi entropies S 2 and S 3 on the 
m = lj6 plateau at low temperatures using the QMG 
extended ensemble method proposed in Ref. 68 (with a 
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small modification where the size of A is dynamically 
varied during the simulation). 

An important first observation is that if one consid¬ 
ers only the smallest (four or five) system sizes, a naive 
fit leads to a negative intercept and thus to an appar¬ 
ent topological phase. One could perhaps interpret this 
unexpected finite-size effect in the following way: close 
enough to the transition to the superfluid phase, the su¬ 
perfluid correlation length can be of the order of, or larger 
than, the system size for small enough samples. Effec¬ 
tively, the system acquires a superfluid component, which 
we expect to contribute to Sq as a'qiA + blog{£A) + d'q. 
This form is derived in Ref. 69 for a system that sponta¬ 
neously breaks a continuous symmetry, in which case the 
logarithm prefactor b is universal (positive and related to 
the number of Goldstone modes), and d'q a non-universal 
constant. Assuming such a form is also valid in a system 
which has effectively both superfluid and crystalline cor¬ 
relations (ie., supersolid), it could result in fits that dis¬ 
play a non-universal, possibly negative (for small enough 
systems, since otherwise b\og{iA) S> |d^| for large enough 
L), constant intercept. Note that it would be interest¬ 
ing to check this behaviour using simpler bosonic models 
that exhibit supersolid behaviour. 

However, when the system is large enough (beyond 
L = 7 in our simulations), we observe a qualitative 
change of behavior in Sq as a function of L. Note that this 
length scale coincides with the one where the stripe struc¬ 
ture factor starts saturating, as observed in Fig. 13. For 
these larger systems, entanglement entropies are com¬ 
patible with a positive intercept dq ~ log(3), although 
the accuracy is not excellent. For instance, fitting our 
S 2 data using (10) and sizes with 7 < L < 10, we get 
dq ~ 0.9, which deviates quite significantly from log(3). 
Indeed, as already advocated, a precise determination of 
dq is rather difficult due to the need of excellent ergodicity 
and hence large autocorrelation times for this quantity. 
Typically, we need more than 10® measurements to get 
reliable data. 


V. CONCLUSION 

In search for magnetization plateaus of spin-liquid na¬ 
ture, we have considered the effect of a magnetic field 
on an anisotropic spin-1/2 model on the Kagome lat¬ 
tice that was introduced in Ref. 10. In the easy-axis 
limit, this model is designed to exhibit several magneti¬ 
zation plateaus, at magnetization per site m = 0,1/6,1/3 
(plus saturation 1/2). The properties of these plateaus 
have been investigated first by mapping the original spin 
model to simpler constrained models on the triangular 
lattice. Focusing on the m = 1/6 plateau, we have in¬ 
vestigated in detail the effective quantum loop model and 
mapped out its phase diagram using a reptation quantum 
Monte-Carlo algorithm. From this, we predict a 3-fold 
degenerate stripe phase, that breaks rotation symmetry, 
to appear on this magnetization plateau. We have also 


observed that the m = 1/6 plateau with the stripe order 
turns into a topological Z 2 spin-liquid plateau when an 
additional repulsion among the dimer segments is added. 

Then, we have performed large-scale numerical simu¬ 
lations using the SSE quantum Monte-Carlo algorithm 
on the original microscopic model. This has confirmed 
the existence of a stripe phase found in the effective 
loop model, using both standard structure-factor mea¬ 
surements and its signature in entanglement entropies. 
In particular, we have observed that, for length scales 
shorter than L ~ 6, the scaling of the block entangle¬ 
ment entropy may be misleadingly interpreted as giving 
a negative intercept dq, while it becomes clearly positive 
for larger systems. This should be used as a caveat in 
other situations. 

When the XY-interaction \Jxy\ increases, there is a 
quantum phase transition to a planar phase (superfluid of 
the equivalent hardcore bosons). While it appears to be 
continuous in our simulations, we believe that it should 
appear weakly first-order on larger lattices, as observed 
for instance in similar physical situations [58-60]. 

We have not investigated in details the m = 1/3 
plateau, where the effective model is given by the stan¬ 
dard quantum dimer model (ie. a single dimer per site) 
on the triangular lattice, that has been widely studied 
in the past [8, 23-29]. In fact, the expected crystal in 
this phase would have a -s/I^ x \/T2 unit cell and an ex¬ 
tremely small order parameter [25] . Civen that numerical 
simulations of the microscopic spin model are more in¬ 
volved (larger systems, many different energy scales etc.), 
it remains challenging to detect such a weak order in the 
direct simulations for the original BFC spin model. 

Before concluding, let us mention possible extensions 
of the model (1). In Sec. Ill, we have seen that there ex¬ 
ists a gapped Z 2 spin liquid phase around the RK point 
of the QLM (see Fig. 7). Therefore, it would be inter¬ 
esting to drive the system toward the spin-liquid phase. 
In fact, it is possible to consider additional interactions 
which would result in a non-zero interaction term V in 
the effective constrained model (5). The four-spin inter¬ 
action Q [Eq. (6)] or the inter-hexagon third-neighbor in- 
teraction J) ' [Eq. (7)], that is not included in the model 
(1), should do the job. This could allow us to reach the 
Z 2 liquid state at both m = 1/3 and m = 1/6 magnetiza¬ 
tion plateaus when ^ 2 Jjy/J 2 , while the transi¬ 

tion between the crystalline and topological phases could 
in principle be monitored by examining the Renyi en¬ 
tropy behavior. This extension of the BFC model would 
certainly be a very interesting playground to investigate 
such topological phases, their detections and their tran¬ 
sitions to conventional symmetry-breaking states. Also 
our construction is not restricted to spin-1/2 systems (or 
to the equivalent hardcore-boson systems) and applies 
to higher-spin systems as well. In fact, we can obtain 
the same effective models (QDM and QLM) for the spin- 
1 version of the model (1) with an additional single-ion 
anisotropy 

As a final remark, let us briefly comment on possible 
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experimental realization of the BFG model. Because of 
the artificial interactions required to build the ice man¬ 
ifold of the BFG model, it is difficult to make direct 
connections between the model we considered and the 
existing Kagome compounds [70]. However, it remains 
an interesting prospect to investigate whether such toy 
models could be realized experimentally in tunable artifi¬ 
cial systems. Indeed, it was suggested very recently that 
the BFG interactions and hexagonal plaquettes could be 
reproduced using a 2D cold ion crystal, with an example 
given for one and two hexagons [71]. 

Note added: While finishing this manuscript, we be¬ 
came aware of the recent preprint by Roychowdhury et 
al. [72] who recently derived the same loop model for 
the equivalent hardcore-boson Hamiltonian, and further 
studied the phase diagram of the loop model with a po¬ 
tential term. While we agree on the phase structure of 
the loop model, the position of the critical point sepa¬ 
rating the stripe from the Z 2 liquid phase is noticeably 
different [{V/t)c —0.3 versus {V/t)c ~ 0.7 ]. 
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Appendix A: Renyi entropy in a discrete symmetry 
breaking phase 

Let us consider g states ipi, each corresponding to one 
of the degenerate states in the thermodynamic limit. 
This would be for instance | titi ■ • ■) and | itit ■ ■ •) in 
an Ising phase, or the three states with down spins occu¬ 
pying one of the three Kagome sublattices (and up spins 
on the other sites, see Fig. 12), in the stripe crystal phase 
discussed in the main text. These g states are always or¬ 
thogonal in the thermodynamic limit. Since they are sim¬ 
ple product states, one can readily show, for any block A 


of the system which complement scales with system size, 
that the reduced density matrix corresponding to a cat 
superposition of these states, IV'cat) = (l/^/s) IV'i): is 
given by: 

= (Al) 

P i 

Note that this equality only holds for the reduced density 
matrix, and not the whole density matrix. 

On any finite system, since these g states are expo¬ 
nentially close in energy, and a QMG simulation would 
sample each I'lpi) with equal probabilities. Hence, one 
would compute a reduced density matrix corresponding 
to a mixed state: 

PQMC = - ^ Pi => PqMC = “ ^ pf (A2) 

^ i ^ i 

which is precisely the same object as Eq. (Al). So our 
QMG simulation is able to access the reduced density 
matrix of the cat state, provided it is ergodic. 

Now, each is excessively simple for each product 
state, as its spectrum contains only one non-zero eigen¬ 
value A = 1 since they are not entangled. Moreover, they 
can be diagonalized simultaneously (they are diagonal 
in the basis in our example), which leads to an en¬ 
tanglement spectrum with non-zero eigenvalues 1 /g with 
degeneracy g. From it, the entanglement entropy (von 
Neumann or Renyi) is simply logg. 

This results looks quite intuitive indeed and was con¬ 
jectured in Ref. [66]. It was already noticed in Ref. [74] 
for a cat state made of g degenerate Ising ground states. 
We emphasize that the same result holds for the mixed 
state obtained in our QMG simulation. For more realistic 
wave-functions, the fluctuations around the ideal product 
state gives and area law, but the constant term will be 
unchanged: 

Sq (A) = -t- log g -f ... (A3) 

As a last remark, we would like to emphasize that in 
other more complex cases (e.g. for states which are not 
one-site product states, such as valence-bond columnar 
states), the discussion is more involved as the constant 
term can now depend explicitly on the form of the cut as 
well as on the Renyi index q. 
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